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We present a fully parameter-free density-functional approach for the accurate description of opti¬ 
cal absorption spectra of insulators, semiconductors and metals. We show that this can be achieved 
within time-dependent current-density-functional theory using a simple dynamical polarization func¬ 
tional. We derive this functional from physical principles that govern optical spectra. Our method 
is truly predictive because not a single parameter is used. In particular, we do not use an ad-hoc 
material-dependent broadening parameter to compare theory to experiment as is usually done. Our 
approach is numerically efficient; the cost equals that of a calculation within the random-phase 
approximation. 

PACS numbers: 71.15.Mb, 71.45Gm, 71.10-w, 71.15Qe 


The calculation of accurate optical absorption spec¬ 
tra within a density-functional approach has been a ma¬ 
jor challenge for almost two decades. While an accu¬ 
rate approach to calculate optical spectra is by solving 
the Bethe Salpeter equation (BSE) [l]-[H, the large nu¬ 
merical effort required precludes its application to large 
systems. Therefore an approach based on the numeri¬ 
cally more affordable time-dependent density-functional 
theory (TDDFT) [fj or time-dependent current-density- 


functional theory (TDCDFT) [6H8| that gives optical 
spectra of similar quality as the BSE is much sought after. 
The price to pay for the gained computational efficiency 
is that the auxiliary Kohn-Sham system of TD(C)DFT is 
governed by an unknown effective potential for which the 
derivation of accurate approximations is highly nontriv¬ 
ial. The main failures of standard approximations such 
as the random-phase approximation (RPA) and the adia¬ 
batic local-density approximation (ALDA) [sj are: 1) the 
underestimation of the onset of the absorption; 2) the un¬ 
derestimation of the intensity of continuum excitons; 3) 
the absence of bound excitons; 4) the absence of Drude 
tails in the spectra of metals. While the first problem 
can be circumvented by replacing the Kohn-Sham eigen¬ 
values by quasiparticle energies obtained within the GW 
approximation 11 (1 111 . the other failures have proven 
to be more complicated to solve. Early attempts to go 
beyond such standard approximations focused on a cor¬ 
rect description of continuum excitons by employing a 
long-range exchange-correlation Dec) kernel either within 
TDCDFT [13 or TDDFT 13. Q. These approaches 
produced the desired result but at the cost of introduc¬ 
ing a material-dependent parameter. The nanoquanta 
kernel 15 does not require such a parameter and 
leads to a correct description of both bound and con¬ 
tinuum excitons. Unfortunately, being derived from the 
BSE, it is almost as expensive to evaluate as solving the 
full BSE. An efficient and parameter-free TDCDFT ap¬ 
proach can be obtained using the Vignale-Kohn func¬ 


tional [20j. While this method describes well the optical 
spectra of metals [21[ it does not describe correctly ex¬ 
citons in semiconductors and insulators 22j. Recently 
several new approaches have been proposed. The boot¬ 
strap method advances an expression for the TDDFT xc 
kernel which has to be obtained from a self-consistent 
procedure [23|]. In Ref. 31 an xc functional is derived 
from a meta-generalized-gradient approximation. Fi¬ 
nally, Ref. 31 proposes an xc functional that is based on 
the jellium-with-gap model. Although these new meth¬ 
ods improve the description of excitonic effects, they have 
several shortcomings: 1) they are not parameter free due 
to an ad-hoc material-dependent broadening parameter 
that is used to compare theory and experiment; 2) they 
are static and therefore do not account for memory ef¬ 
fects precluding a parameter-free description of the finite 
width of Drude tails and bound excitons; 3) they require 
the calculation of the full Kohn-Sham density-density re¬ 
sponse function, an 0(N 4 ) calculation (with N the num¬ 
ber of atoms or electrons in the unit cell). 


In this work we present an efficient 0(N 3 ) density func¬ 
tional approach that does not suffer from the shortcom¬ 
ings mentioned above. It is parameter free and dynamic, 
and it produces accurate optical spectra for metals, semi¬ 
conductors, and insulators. The method we propose is 
both simple and elegant and is characterized by the fol¬ 
lowing two steps: 1) We calculate the macroscopic di¬ 
electric function using an efficient TDCDFT approach 
31 [ 27 } within the RPA (including a scissors shift ob¬ 
tained from GW ); 2) We apply a simple dynamical po¬ 
larization functional (see Eq. (fT4l) below) that accounts 
for continuum and bound excitons as well as Drude tails. 
We will now give all the details of our approach. We use 
Hartree atomic units throughout. 


The optical response of a solid can be obtained as the 
response to a homogeneous electric field, i.e., a field with 
vanishing wave vector q. However, for <f = 0, the periodic 
density is not enough to uniquely describe the response; 
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one also needs the macroscopic polarization [281 ]. This is 
an immediate consequence of the use of periodic bound¬ 
ary conditions. Macroscopic effects due to the surface 
density have to be accounted for via the macroscopic po¬ 
larization. An elegant and efficient approach is to use the 
current density as the fundamental quantity of a DFT 
for time-dependent phenomena. The self-consistent per¬ 
turbing potential is a unique functional of the periodic 
current density since the latter contains both the peri¬ 
odic density of the bulk (via the continuity relation) and 
the macroscopic polarization, which is obtained as a bulk 
property using the definition, 


M = ^ Jydr^ir^). 


(1) 


Here V is the volume of the unit cell. The limit q —» 0 
can thus be performed explicitly, by setting q = 0, and 
hence knowledge of the current density suffices and no 
explicit calculation of Kohn-Sham response functions is 
needed. The macroscopic polarization is induced by a 
macroscopic electric field E mac (ui) which comprises both 
the externally applied field and the macroscopic induced 
electric field. The constant of proportionality is the elec¬ 
tric susceptibility tensor Xe{w) which is defined as 


Pmac {^0 — Xe(^) ' E rnac {uj'). (2) 

The macroscopic dielectric tensor ^m(w) can be obtained 
from the electric susceptibility according to 

^m(w) = 1 + 47TXe(w). (3) 


Therefore, for a given E mac (oj ), knowledge of the induced 
current suffices to calculate 'c'm(w) from Eqs. ©-©. 

Within Kohn-Sham linear response the current density 
is given by 


kernel f xc {f, f, w): 



Substitution of Eq. © into Eq. ©, and using Emac(w) = 
iuAmac{oj), reveals that P m ac(u) is linear in the macro¬ 
scopic Kohn-Sham electric field, E mac (u})+E%£ ac (w). Fol¬ 
lowing Ref. jlH , we define the auxiliary susceptibility x° 
according to 

Pmacico) = t(w) • \E mac {u) + E% ac (u)\. (6) 

where x° i s the susceptibility obtained from a calcula¬ 
tion with E^ ac (uj) = 0. We note that the contribu¬ 
tions of 5v*?c{r, oj) to x° are fully accounted for. Since 
S u^f c c (r, w) is itself a functional of Sj(f,uj) this is done 
within a self-consistent field (SCF) calculation. We will 
now show that the macroscopic xc effects can be ac¬ 
counted for post-SCF. Comparison of Eqs. © and © 
leads to the following relation between Xe and Xe, 

(RTV) - [XePV)) ■ Pma C {u) = E % a ». (7) 

If one chooses Ef^ ac {oj) equal to zero as is usually done 
we obtain Xe(w) = Xe( w )- We can now go beyond this 
approximation and obtain a polarization functional for 
£ , ^ l c ac (w) by neglecting microscopic current components 
in Eq. ©,*.e., we replace Sj(f. oj) by its unit-cell average. 
We obtain 

^mac(^) = &(<*>) ' Pmaci^) (8) 

where we used Eq. © and in which we defined 

a(w) = -y J dr J dr' f xc {r,r*,w). (9) 

Substitution of Eq. © into Eq. © leads to 


5j(r,u) = J dr ,i jt s ’ n {r,r , ,u)[A rnac (uj) + $£ ac (? ,u)} 

+ J d?x s3p {?y , u), (4) 

where A mac (oj ) = E mac (ui)/(iu). We note that the 
Kohn-Sham current-current and current-density response 
functions (x s ’-”(w) and x s,3/, (o;), respectively) are never 
needed explicitly since we use a sum-over-states repre¬ 
sentation. Therefore the r' integrals can be evaluated 
independently of r. In Eq. © we used the micro¬ 
scopic Coulomb gauge in which the microscopic poten¬ 
tial 6v^fc(r, oj) is lattice periodic and contains the peri¬ 
odic part of the Hartree and longitudinal xc contributions 
[26j . All remaining xc contributions are included in the 
the vector potential A xc (r,oj). Here we neglect micro¬ 
scopic transverse xc contributions and therefore only its 
macroscopic part, A^ ac (oj), enters Eq. ©. It is related 
to the induced current through the TDCDFT tensor xc 


Ke]- 1 M = K°]- 1 H-S( w ) (10) 

Therefore, for a given a(oj), we can simply calculate 
X e (w) from x°M- 

Here we will use the RPA (Sv^f^ = Sv^ ic ) to calcu¬ 
late Xe(w), i.e., Xe = X PPA - Since continuum excitons 
are underestimated and bound excitons and Drude tails 
are absent in RPA optical spectra, we will include these 
effects through a(oj). We first derive a static a from 
the physical principles that govern strongly bound ex¬ 
citons in solids. To simplify our argumentation we will 
assume for the moment that x e (w), X PPA ( UJ ) and a(w) 
are isotropic, i.e., a.ij(oj) = a(ui)5ij, etc.. Following ar¬ 
guments similar to those given in Ref. [29], we derive in 
the supplemental material [13] that the exact constraint 
to have a bound exciton is that for some oj = ojbe below 
the band gap the following relation holds for a(oj), 

Re[a(w be )] = l/x PPA (vbe) and Im[a(u;6 e )] ~ 0. 

( 11 ) 



However, the application of the above expression requires 
the knowledge of u>b e which is unknown. Therefore, we 
would like to rewrite the above expression in terms of 
known quantities. In the supplemental material [30j we 
derive such an expression for the case of a wide-gap in¬ 
sulator whose optical spectra is dominated by a bound 
exciton. The result is a static a given by 


1 

Xe PA (= 0)£m PA (w = 0)' 


( 12 ) 


Although the derivation applies to wide-gap insulators, 
Eq. (fT^l) has the additional advantage that a is propor¬ 
tional to [eA/’ A (0)] _1 . It has previously been shown nu¬ 
merically that such a proportionality leads to good op¬ 
tical spectra for semiconductors [lj, |31[. The physical 
reason is that e^ p " 4 (0) is a measure of the amount of 
screening of the electron-hole interaction. Screening ef¬ 
fects are more important in semiconductors than in wide- 
gap insulators. We note that Eq. m has a similar form 
as the bootstrap kernel [23] but that no self-consistent 
procedure is needed to calculate it. Unfortunately a in 
Eq. (fl2l) is static and will therefore not be able to account 
for the finite width of bound excitons and Drude tails. 
For this reason we add to Eq. (fl2l) Yvk(u), the long- 
range part of the dynamical Vignale-Kohn functional [201 ] 
in which we replace the current density by its unit-cell 

average [ 13 ]: 


Y VKM =lfdr y ^ 


+ 


v V 

Vpo(r) ® Vpo(r) 


(Vp 0 (r) • Vpo(r) r ^ 
Of .-A JxcTyPi 


AA 

ui) I 


fxcL(p,v) ~ fxcT(p,v) - 


dp 


a2 


(13) 


Here fxcL^T)^) is the longitudinal (transverse) xc ker¬ 
nel of the homogeneous electron gas, e xc is the xc en¬ 
ergy per volume of the homogeneous electron gas, po(r) 
is the ground-state density and p is its average in the unit 
cell. For f xc nr )( w ) we use the parametrization of Ref. 
[32 ] in the QVA approximation 211. The Vignale-Kohn 
functional is the exact functional of a slightly inhomo¬ 
geneous electron gas and describes correctly the optical 
spectra of metals 2lj. For this reason Eq. (TT3l) is com¬ 
plementary to Eq. m which tends to zero for metal¬ 
lic systems because screening is complete and therefore 
[ e M PA ( w = 0)] -1 0- Yvk(v) will account for the 

Drude tails and finite width of bound excitons which 
are exactly the features that are absent in Eq. (1121) . 
Moreover, as we show in the supplemental material [30(, 

AA 

Yvk{ w) will not much influence the spectra of semicon¬ 
ductors. We finally obtain the following approximation 
for a(uj), 


AA / \ _ [AA 

a{cu) = [e 


++RPA ( 


■\-l\ttRPA 


M 


®r A m~ l +YvK m (i4) 




FIG. 1: (Color online) The optical absorption spectra of bulk 
silicon and GaP. Solid line (black): polarization functional 
(PF); Dashed line (red): RPA; Dotted line (blue): experiment 
from Ref. fH (Si) and Ref. gjj (GaP). 


where we generalized Eq. m to a tensor form. Since 
[^(O)]- 1 and [xf PA (0)] _1 commute, the order of the 
multiplication is irrelevant. We note that Eq. satis¬ 
fies the Kramers-Kronig relations. Equation (1141) is the 
main result of this work. 

We will now demonstrate our approach by applying it 
to the calculation of the optical spectra of several mate¬ 
rials. We briefly outline the full procedure. The ground- 
state calculations are done within the local-density ap¬ 
proximation (LDA) [33| and we use LDA lattice param¬ 
eters in all calculations. We apply a scissors operator to 
shift the unoccupied bands and we modify the current 
operator accordingly to guarantee that exact constraints 
such as the continuity equation remain satisfied. The 
energy shift is calculated with the GW method 0,0 
and is equal to the Go Wo correction for the direct band 
gap at the T point. We calculated this correction us¬ 
ing the effective-energy technique [34]-[36|. which leads 
to a speed-up of an order of magnitude. The calcula¬ 
tions were done with the Abinit code [37]. We imple¬ 
mented the polarization functional in a modified version 
of the Amsterdam Density Functional (ADF) code [38- 
0 - We use the TZ2P (triple-^ + 2 polarization func¬ 
tions) basis set provided by ADF. The fc-space integrals 
are done analytically using a Lehmann-Taut tetrahedron 
scheme [4lJ - Since we do not include effects due to 
electron-phonon coupling the spectra obtained with the 
above approach are predictions of the optical spectra at 
low temperature where electron-electron scattering dom¬ 
inates electron-phonon scattering. For this reason we will 
compare our calculated spectra with spectra measured at 
low temperature where available. 

In Fig. [Qwe report the optical absorption spectra of 
bulk silicon and bulk GaP obtained with our polariza¬ 
tion functional and compare it to the RPA spectra and 
to experimental results obtained at low temperature (15 
Kelvin). Silicon and GaP are typical examples of mate- 































i 1 ~ 

* 1 1 1 T 1 ^ 

LiF - 

- 

- 

- 

- 

- 

- 

: 

- 


■L 


0 I 1 — , _C 

10 12 14 16 18 

CO(eV) 




FIG. 2: The optical absorption spectra of solid Argon and 
LiF. Solid line (black): polarization functional (PF); Dashed 
line (red): RPA; Dotted line (blue): experiment from Ref. 
0 (Ar) and Ref. [H (LiF). 


FIG. 3: The optical absorption spectra of diamond and cop¬ 
per. Solid line (black): polarization functional (PF); Dashed 
line (red): RPA; Dotted line (blue): experiment from Ref. 
(4til (Diamond) and Refs. [47j and |48j] (Cu). 


rials for which the RPA strongly underestimates the first 
peak which appears in the experimental spectra around 
3.4 eV (Si) and 3.8 eV (GaP). Our polarization functional 
solves this problem by including the necessary excitonic 
effects and the first peak compares well with experiment 
both in position and magnitude. Overall, the spectra are 
very close to experiment with the exception of the peak 
around 5.2 eV in the spectrum of GaP which is overesti¬ 
mated. 

In Fig.[2]we show the optical absorption spectra of solid 
argon and LiF obtained with our polarization functional 
and compare it to the RPA spectra and to experimental 
results. Solid argon and LiF are typical materials that 
exhibit strongly bound excitons. We see that these ex- 
citons which appear in the experimental spectra around 
12 eV (Ar) and 12.5 eV (LiF) are completely absent in 
the RPA spectra. Our polarization functional describes 
these bound excitons and also accurately reproduces their 
position. The magnitude of the peaks is overestimated 
with respect to experiment. This discrepancy is prob¬ 
ably due to electron-phonon broadening in the experi¬ 
ments and to the fact that density-functional approaches 
tend to overestimate these peaks [23|. ;25|. Finally, we 
verify a posteriori that a given in Eq. m is indeed a 
good approximation to Eq. m for wide-gap insulators 
exhibiting a dominant bound exciton. With Eq. m we 
obtain auF = 8.9 and QAr = 11-8, while Eq. m gives 
QTiF = 9.2 and «A r = 12.2. 

In Fig. [3] we report the optical absorption spectra of 
diamond and copper obtained with the polarization func¬ 
tional and compare it to the RPA spectra and to experi¬ 
mental results. Diamond is another typical test case since 
the RPA spectrum is quite different from the experimen¬ 
tal spectrum. Due to the absence of excitonic effects the 
RPA spectrum has too much weight at high energy. With 
our polarization functional the spectral weight is shifted 
to lower energy and we obtain a very good agreement 


with experiment. While the RPA spectrum of copper 
accurately reproduces the part of the spectrum which is 
due to interband transitions, the Drude tail at low en¬ 
ergy, which is due to intraband transitions, is completely 
absent. Our polarization functional accurately describes 
the Drude tail while maintaining the good agreement for 
the inter band part. 

In conclusion, we presented the first fully parameter- 
free density-functional approach that gives accurate op¬ 
tical spectra for insulators, semiconductors and metals 
alike. Our approach is therefore truly predictive and due 
to its numerical efficiency opens the way for the predic¬ 
tion of optical spectra of large systems. 
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